R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #################################################################
> # "Democracy, Inequality and Antitrust"                         #
> #   By Michael O Allen, Kenneth Scheve, and David Stasavage     #
> #                                                               #
> #   Last modified: 12/21/2024                                   #
> #################################################################
> 
> library(data.table)
data.table 1.15.0 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********
> library(ggplot2)
Keep up to date with changes at https://tidyverse.org/blog/
> library(cowplot)
> library(PanelMatch)
> library(fixest)
> library(marginaleffects)
> library(interflex)
## Syntax has changed since v.1.2.1.

## See http://bit.ly/interflex for more info.
## Comments and suggestions -> zyliu2020@uchicago.edu.

> library(xtable)
> library(modelsummary)
> 
> library(here)
here() starts at /Users/moda/Desktop/democracy inequality antitrust jop desk
> 
> 
> # Function ----------------------------------------------------------------
> # modelsummary not printing adjusted R2 for some reason... so need these:
> getGroups <- Vectorize(function (x) x$fixef_sizes[['code_cowc']])
> 
> # Adding commas to things
> f <- function (x) format(x, big.mark = ',')
> 
> 
> # Table 1 -----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> ### Analyses ------------------------------------------------------------
> 
> t1_1 <- feols(p3_CLI_DV ~  dem_lag1 | code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 210 observations removed because of NA values (LHS: 25, RHS: 210, Fixed-effects: 25).
> 
> t1_2 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 895 observations removed because of NA values (LHS: 25, RHS: 895, Fixed-effects: 25).
> 
> t1_3 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 901 observations removed because of NA values (LHS: 25, RHS: 901, Fixed-effects: 25).
> 
> t1_4 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 915 observations removed because of NA values (LHS: 25, RHS: 915, Fixed-effects: 25).
> 
> t1_5 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 +
+               crisis_3_percent_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 916 observations removed because of NA values (LHS: 25, RHS: 916, Fixed-effects: 25).
> 
> 
> ### Gen table -------------------------------------------------------------
> 
> coefNames <- c("dem_lag1" = "Democracy",
+                "top1Share_interp_lag1_ln" = "ln Top 1% Income Share",
+                "dem_lag1:top1Share_interp_lag1_ln" = "ln Top 1% Income Share x Democracy",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list(t1_1, t1_2, t1_3, t1_4, t1_5)
> names(models) <- paste0("(", 1:length(models), ")")
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", 5),
+                                 "Country & Period FE?", rep("Yes", 5),
+                                 "Countries", getGroups(models)),
+                               ncol = 6, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country",
+          output = here('output', 'tables', 'table-1.tex')
+ )
Warning message:
To compile a LaTeX document with this table, the following commands must be placed in the document preamble:

\usepackage{booktabs}
\usepackage{siunitx}
\newcolumntype{d}{S[
    input-open-uncertainty=,
    input-close-uncertainty=,
    parse-numbers = false,
    table-align-text-pre=false,
    table-align-text-post=false
 ]}

To disable `siunitx` and prevent `modelsummary` from wrapping numeric entries in `\num{}`, call:

options("modelsummary_format_numeric_latex" = "plain")
 This warning appears once per session. 
> 
> 
> ## Scaled estimate ----------------------------------------------------------------------
> 
> m_scaled_dv <- feols(scale(p3_CLI_DV) ~  dem_lag1 * top1Share_interp_lag1_ln +
+                        gdp_lag1_ln +
+                        gdppc_lag1_ln +
+                        trade_openness_lag1 +
+                        crisis_3_percent_lag1 | 
+                        code_cowc + p3,
+                      data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 916 observations removed because of NA values (LHS: 25, RHS: 916, Fixed-effects: 25).
> 
> # Log-point increase in inequality reduces association 
> #   between democracy and CLI by .674 SDs
> m_scaled_dv$coefficients['dem_lag1:top1Share_interp_lag1_ln']
dem_lag1:top1Share_interp_lag1_ln 
                       -0.6740273 
> 
> # Figure 2 ----------------------------------------------------------------
> 
> # Uses estimates from Table 1 above
> cme_data <- plot_cme(t1_5, 
+                      variables = "dem_lag1", 
+                      condition = "top1Share_interp_lag1_ln",
+                      draw = FALSE)
> ggplot(cme_data,
+        aes(x = top1Share_interp_lag1_ln, y = estimate)) +
+   geom_hline(yintercept = 0, color = 'red', size = .5, linetype = 'dashed') +
+   geom_line() +
+   geom_ribbon(aes(ymin = conf.low, 
+                   ymax = conf.high), alpha = .2) +
+   labs(x = 'ln Top 1% Income Share',
+        y = 'Marginal Effect of Democracy on CLI') +
+   scale_x_continuous(limits = c(1.5, 3.5)) +
+   scale_y_continuous(limits = c(-.25,.5),
+                      breaks = seq(-.2, .5, .1)) +
+   theme_minimal(base_size = 14) +
+   theme(legend.position = 'none',
+         panel.grid.minor.x = element_blank(),
+         panel.grid.minor.y = element_blank())
Warning messages:
1: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated. 
2: Removed 20 rows containing missing values or values outside the scale range (`geom_line()`). 
> 
> ggsave(here('output', 'figures/figure-2.png'),
+        width = 6, height = 4, units = 'in', dpi = 450)
Warning message:
Removed 20 rows containing missing values or values outside the scale range (`geom_line()`). 
> 
> 
> # Figure 3 ----------------------------------------------------------------
> 
> # Load the function to plot
> source(here('scripts', 'plot-function.R'))
> 
> ## Data -------------------------------------------------------------------
> 
> df_f3 <- fread(here('data', 'cli_yearly.csv'))
> 
> df_f3[ , cid := as.integer(factor(code_cowc))]
> df_f3[ , gdp_ln := log(gdp)]
> df_f3[ , gdppc_ln := log(gdppc)]
> 
> vars <- c('cli_overall_norm',
+           'gdp_ln',
+           'gdppc_ln',
+           'trade_openness',
+           'top1Share_interp_ln')
> 
> df_f3[ , paste0(vars, '_lm3') := (frollmean(.SD, 3, align = 'right',
+                                                na.rm = T, hasNA = T)),
+           .SDcols = vars,
+           by = code_cowc]
> 
> ## Create subsets ----------------------------------------------------------
> # Low inequality transitions only
> # Drop transitions + decade that occurred with above median levels of inequality or unknown
> 
> df_f3$lowInequality_subset_drop <- 1
> 
> for (r in which((df_f3$top1Share_interp >= median(df_f3$top1Share_interp, na.rm = T) |
+                  is.na(df_f3$top1Share_interp)) 
+                 & df_f3$dem_trans == 1)) {
+   df_f3[ r:(r+10), lowInequality_subset_drop := 0]
+ }
> 
> 
> # High inequality subset
> # Drop transitions that occurred with below median levels of inequality or unknown
> 
> df_f3$highInequality_subset_drop <- 1
> 
> for (r in which((df_f3$top1Share_interp < median(df_f3$top1Share_interp, na.rm = T) |
+                  is.na(df_f3$top1Share_interp)) 
+                 & df_f3$dem_trans == 1)) {
+   df_f3[ r:(r+10), highInequality_subset_drop := 0]
+ }
> 
> # Analysis -----------------------------------------------------------------
> 
> ## Low inequality -------------------------------------------------------------------
> 
> data_polity_low <- as.data.frame(df_f3[lowInequality_subset_drop == 1 & year >= 1945])
> 
> pm_lowInequality <- PanelMatch(lag = 5, lead = 0:10, 
+                                time.id = 'year', unit.id = 'cid', 
+                                treatment = 'dem_bin', 
+                                outcome.var = 'cli_overall_norm',
+                                refinement.method = 'CBPS.weight',
+                                covs.formula = ~ 
+                                  lag(cli_overall_norm_lm3, 1) +
+                                  lag(gdp_ln_lm3, 1) +
+                                  lag(gdppc_ln_lm3, 1) +
+                                  lag(trade_openness_lm3, 1) +
+                                  lag(top1Share_interp_ln_lm3, 1),
+                                qoi = 'att', 
+                                match.missing = FALSE,
+                                forbid.treatment.reversal = FALSE,
+                                data = data_polity_low)
Warning message:
In panel_match(lag, time.id, unit.id, treatment, refinement.method,  :
  non-numeric data exists. Only numeric (including binary) data can be used for refinement and calculations
> 
> set.seed(1)
> pe_low <- PanelEstimate(pm_lowInequality, number.iterations = 2000,
+                         data = data_polity_low)
> low_est <- round(summary(pe_low)$summary[, 1:2] ,3)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
> 
> ## High inequality ------------------------------------------------------------------
> 
> data_polity_high <- as.data.frame(df_f3[highInequality_subset_drop == 1 & year >= 1945])
> 
> pm_highInequality <- PanelMatch(lag = 5, lead = 0:10, 
+                                 time.id = 'year', unit.id = 'cid', 
+                                 treatment = 'dem_bin', 
+                                 outcome.var = 'cli_overall_norm',
+                                 refinement.method = 'CBPS.weight',
+                                 covs.formula = ~ 
+                                   lag(cli_overall_norm_lm3, 1) +
+                                   lag(gdp_ln_lm3, 1) +
+                                   lag(gdppc_ln_lm3, 1) +
+                                   lag(trade_openness_lm3, 1) +
+                                   lag(top1Share_interp_ln_lm3, 1),
+                                 qoi = 'att', 
+                                 match.missing = FALSE, 
+                                 forbid.treatment.reversal = FALSE,
+                                 data = data_polity_high)
Warning message:
In panel_match(lag, time.id, unit.id, treatment, refinement.method,  :
  non-numeric data exists. Only numeric (including binary) data can be used for refinement and calculations
> 
> # Used to make Table C.1
> set.seed(1)
> pe_high <- PanelEstimate(pm_highInequality, number.iterations = 2000,
+                          data = data_polity_high)
> high_est <- round(summary(pe_high)$summary[, 1:2] ,3)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
> 
> ## Create plots ------------------------------------------------------------
> 
> nicePMplot(pm_lowInequality, data_polity_low,
+            yupper = .4, ylower = -0.21,
+            seed = 1)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
$plot

$estimates
       estimate        2.5%     97.5%          5%       95%
t+0  0.05964276 -0.02417865 0.2212040 -0.01977147 0.1854654
t+1  0.09498239 -0.03913495 0.2488507 -0.02137613 0.2138525
t+2  0.08788012 -0.03379333 0.2295329 -0.01799769 0.1976716
t+3  0.19928400  0.03778300 0.3822526  0.05413377 0.3463673
t+4  0.18320736  0.01617541 0.3614566  0.03436793 0.3278199
t+5  0.17834001  0.01435881 0.3542017  0.02919702 0.3185453
t+6  0.19283228  0.01925463 0.3606922  0.04707318 0.3276121
t+7  0.19899114  0.02309248 0.3584873  0.05553313 0.3268149
t+8  0.19810213 -0.01110577 0.3689941  0.03173036 0.3370448
t+9  0.20905419  0.01624491 0.3793849  0.05445201 0.3471979
t+10 0.19931291  0.01209385 0.3673414  0.04588501 0.3354596

> 
> ggsave(here('output', 'figures', 'figure-3a.png'),
+        width = 6, height = 4, units = 'in', dpi = 450)
> 
> nicePMplot(pm_highInequality, data_polity_high,
+            yupper = .4, ylower = -0.21,
+            seed = 1)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
$plot

$estimates
         estimate        2.5%        97.5%          5%          95%
t+0  -0.009088547 -0.02683497 -0.002814669 -0.02088075 -0.003350689
t+1   0.034545932 -0.04438958  0.163132556 -0.03275127  0.132079257
t+2   0.023816911 -0.06884806  0.149102257 -0.05037233  0.121260566
t+3   0.067094577 -0.05683389  0.216316367 -0.03932350  0.184272010
t+4   0.049443974 -0.09347513  0.194721781 -0.06516368  0.165288396
t+5   0.045801239 -0.10181139  0.189609399 -0.07227939  0.162079838
t+6   0.066126299 -0.11646630  0.230938461 -0.08009329  0.209464900
t+7   0.069804076 -0.12580096  0.220333782 -0.08391908  0.199352920
t+8   0.056237963 -0.15287574  0.206523209 -0.10625709  0.183057285
t+9   0.041432974 -0.17387139  0.187224401 -0.12418238  0.161440462
t+10  0.025232022 -0.20578685  0.174651267 -0.15096446  0.147361703

> 
> ggsave(here('output', 'figures', 'figure-3b.png'),
+        width = 6, height = 4, units = 'in', dpi = 450)
> 
> # Table C1 ----------------------------------------------------------------
> 
> # This formats the estimates and SEs from PanelMatch into a table.
> 
> low_est <- t(matrix(c(low_est[ , 1], paste0( "(", low_est[, 2], ")")),
+        ncol = 2))
> low_est <- matrix(low_est, byrow = TRUE)
> 
> high_est <- t(matrix(c(high_est[ , 1], paste0( "(", high_est[, 2], ")")),
+                     ncol = 2))
> high_est <- matrix(high_est, byrow = TRUE)
> 
> t <- rep("", 22)
> t[seq(1, 21, 2)] <- 0:10
> 
> tc1 <- xtable(cbind(t, low_est, high_est))
> 
> print(tc1,
+       include.rownames = FALSE,
+       file = here('output', 'tables', 'table-c1.tex'))
> # Table A1 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> 
> ## Analyses -------------------------------------------------------------
> 
> 
> ta1_1 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 911 observations removed because of NA values (LHS: 25, RHS: 911, Fixed-effects: 25).
> 
> ta1_2 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 917 observations removed because of NA values (LHS: 25, RHS: 917, Fixed-effects: 25).
> 
> ta1_3 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 930 observations removed because of NA values (LHS: 25, RHS: 930, Fixed-effects: 25).
> 
> ta1_4 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 +
+                  crisis_lag1  
+                  | code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 1,252 observations removed because of NA values (LHS: 25, RHS: 1,252, Fixed-effects: 25).
> 
> ta1_5 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 +
+               crisis_3_percent_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 931 observations removed because of NA values (LHS: 25, RHS: 931, Fixed-effects: 25).
> 
> 
> ## Generate Table A1 ------------------------------------------------------
> 
> coefNames <- c("top1Share_lag1_ln" = "ln Top 1\\% Income Share\\textsubscript{Observed}",
+                "dem_lag1" = "Democracy",
+                "dem_lag1:top1Share_lag1_ln" = "ln Top 1\\% Income Share\\textsubscript{Observed} $\\times$ Democracy",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_lag1" = "Economic Crisis\\textsubscript{RR}",
+                "crisis_3_percent_lag1" = "Economic Crisis\\textsubscript{GDP -3\\%}"
+ )
> 
> models <- list("(1)" = ta1_1, 
+                "(2)" = ta1_2, 
+                "(3)" = ta1_3, 
+                "(4)" = ta1_4,
+                "(5)" = ta1_5)
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = 6, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          escape = FALSE,
+          title = "Re-run of Table \\ref{tab:crossnat_main} with no interpolated inequality data",
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          output = here('output', 'tables', 'table-a1.tex')
+ )
> 
> 
> # Table A2 ----------------------------------------------------------------
> 
> ta2 <- df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ]
> ta2 <- ta2[unlist(t1_5$obs_selection), ]
> setDT(ta2)
> 
> ta2 <- ta2[ , .(obs = .N,
+                   firstPeriod = p3[1]) , by = country]
> 
> periods <- names(attributes(df_t1$p3)$labels)
> 
> ta2$firstYear <- substring(periods[ta2$firstPeriod], 2,5)
> 
> ta2 <- xtable(ta2[ , c('country', 'obs', 'firstYear')])
> print(ta2,
+       file = here('output', 'tables', 'table-a2.tex'))
> 
> rm(ta2, periods)
> 
> # Figure B1 ---------------------------------------------------------------
> 
> ## Function ---------------------------------------------------------------
> 
> # This function collapses the yearly data into different time periods
> #   and estimates the models for Figure B1.
> 
> altDuration_CME <- function(df) {
+   results <- data.table()
+   
+   for (duration in 1:8) {
+     
+     df_temp <- df
+     
+     if (duration != 7) {
+       df_temp$period <- cut(df_temp$year, seq(1890, 2010, by=duration),
+                             right=T)  
+     } else  { # To ensure we include 2010 in the analysis
+       df_temp$period <- cut(df_temp$year, seq(1891, 2010, by=duration),
+                             right=T)
+     }
+     
+     df_temp <- df_temp[!is.na(period)]
+     
+     setorder(df_temp, code_cowc, year)
+     
+     df_temp[ , DV := cli_overall_norm[1] , by = .(code_cowc, period)]
+     
+     df_temp <- df_temp[, .(DV = unique(DV),
+                            top1Share_interp = mean(top1Share_interp, na.rm = TRUE),
+                            dem = mean(dem_bin, na.rm = TRUE),
+                            gdppc = mean(gdppc, na.rm = TRUE),
+                            gdp = mean(gdp, na.rm = TRUE),
+                            crisis_3_percent = mean(crisis_3_percent, na.rm = TRUE),
+                            trade_openness = mean(trade_openness, na.rm = TRUE),
+                            firstYear = min(year),
+                            dem_trans = sum(dem_trans, na.rm = TRUE),
+                            top1_n = unique(top1_n)),
+                        by = .(code_cowc, period)]
+     
+     setorder(df_temp, code_cowc, period)
+     
+     df_temp[ , c('top1Share_lag1',
+                  'dem_lag1',
+                  'gdppc_lag1',
+                  'gdp_lag1',
+                  'crisis_3_percent_lag1',
+                  'trade_openness_lag1',
+                  'dem_trans_lag1') := shift(.SD, 1, NA, 'lag'),
+              .SDcols = c('top1Share_interp', 'dem', 'gdppc', 'gdp', 
+                          'crisis_3_percent', 'trade_openness',
+                          'dem_trans'),
+              by = code_cowc]
+     
+     df_temp[ , gdp_lag1 := log(gdp_lag1)]
+     df_temp[ , gdppc_lag1 := log(gdppc_lag1)]
+     df_temp[ , top1Share_lag1 := log(top1Share_lag1*100)]
+     df_temp[ , p_int := as.integer(period)]
+     
+     df_temp <- df_temp[top1_n >= 20]
+     
+     model <- feols(DV ~ top1Share_lag1 * dem_lag1 
+                    + gdppc_lag1 
+                    + gdp_lag1 
+                    + trade_openness_lag1
+                    + crisis_3_percent_lag1
+                    | code_cowc + period,
+                    data = df_temp)
+     
+     cme_data <- plot_cme(model, 
+                          variables = "dem_lag1", 
+                          condition = "top1Share_lag1",
+                          draw = FALSE)
+     
+     
+     cme_data$period_length <- duration
+     results <- rbind(results, cme_data[, c("period_length", "estimate", 
+                                            "conf.low", "conf.high",
+                                            "top1Share_lag1")])
+   }
+   
+   return(results)
+ }
> 
> 
> ## Load data --------------------------------------------------------------
> 
> df_yearly <- fread(here('data', 'cli_yearly.csv'))
> 
> ## Generate Figure B1 -----------------------------------------------------
> 
> results <- altDuration_CME(df_yearly)
NOTE: 4,514 observations removed because of NA values (RHS: 4,514).
NOTE: 2,287 observations removed because of NA values (RHS: 2,287).
NOTE: 1,580 observations removed because of NA values (RHS: 1,580).
NOTE: 1,174 observations removed because of NA values (RHS: 1,174).
NOTE: 948 observations removed because of NA values (RHS: 948).
NOTE: 822 observations removed because of NA values (RHS: 822).
NOTE: 741 observations removed because of NA values (RHS: 741).
NOTE: 649 observations removed because of NA values (RHS: 649).
> 
> results$period_length <- factor(results$period_length, 
+                                 levels = 1:8,
+                                 labels = paste(1:8, 'year periods'))
> 
> ggplot(results,
+        aes(x = top1Share_lag1, y = estimate)) +
+   geom_line() +
+   geom_hline(yintercept = 0, color = 'red', size = .5, linetype = 'dashed') +
+   geom_ribbon(aes(ymin = conf.low, 
+                   ymax = conf.high), alpha = .2) +
+   labs(x = 'ln Top 1% Income Share',
+        y = 'Conditional Marginal Effect of Democracy on CLI') +
+   scale_x_continuous(limits = c(1.5, 3.5)) +
+   theme_minimal(base_size = 14) +
+   facet_wrap(~period_length) +
+   theme(legend.position = 'none')
Warning message:
Removed 20 rows containing missing values or values outside the scale range (`geom_line()`). 
> 
> ggsave(here('output', 'figures', 'figure-b1.pdf'),
+        width = 7, height = 7)
Warning message:
Removed 20 rows containing missing values or values outside the scale range (`geom_line()`). 
> 
> rm(altDuration_CME, results)
>
>
> # Table B1 ----------------------------------------------------------------
> 
> df_yearly <- fread(here('data', 'cli_yearly.csv'))
> 
> df_yearly[ , gdp_ln := log(gdp)]
> df_yearly[ , gdppc_ln := log(gdppc)]
> 
> tb1_1 <- feols(f(cli_overall_norm) ~ dem_bin * top1Share_interp_ln 
+                | code_cowc + year,
+                panel.id = ~code_cowc + year,
+                df_yearly[ year >= 1985])
NOTE: 2,537 observations removed because of NA values (LHS: 198, RHS: 2,430).
> 
> tb1_2 <- feols(f(cli_overall_norm) ~ dem_bin * top1Share_interp_ln 
+                + gdp_ln
+                + gdppc_ln
+                | code_cowc + year,
+                panel.id = ~code_cowc + year,
+                df_yearly[ year >= 1985])
NOTE: 2,576 observations removed because of NA values (LHS: 198, RHS: 2,472).
> 
> tb1_3 <- feols(f(cli_overall_norm) ~ dem_bin * top1Share_interp_ln 
+                + gdp_ln
+                + gdppc_ln
+                + trade_openness
+                | code_cowc + year,
+                panel.id = ~code_cowc + year,
+                df_yearly[ year >= 1985])
NOTE: 2,581 observations removed because of NA values (LHS: 198, RHS: 2,477).
> 
> tb1_4 <- feols(f(cli_overall_norm) ~ dem_bin * top1Share_interp_ln 
+                + gdp_ln
+                + gdppc_ln
+                + trade_openness
+                + crisis_3_percent
+                | code_cowc + year,
+                panel.id = ~code_cowc + year,
+                df_yearly[ year >= 1985])
NOTE: 2,581 observations removed because of NA values (LHS: 198, RHS: 2,477).
> 
> coefNames <- c("dem_bin" = "Democracy",
+                "top1Share_interp_ln" = "ln Top 1\\% Income Share",
+                "dem_bin:top1Share_interp_ln" = "ln Top 1\\% Income Share $\\times$ Democracy",
+                "gdp_ln" = "ln GDP",
+                "gdppc_ln" = "ln GDP per cap.",
+                "trade_openness" = "Trade Openness",
+                "crisis_3_percentTRUE" = "Economic Crisis"
+ )
> 
> models <- list("(1)" = tb1_1, 
+                "(2)" = tb1_2, 
+                "(3)" = tb1_3, 
+                "(4)" = tb1_4)
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("1 year", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-b1.tex')
+ )
> 
> # Table B2 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> 
> ## Analyses ---------------------------------------------------------------
> 
> tb2_1 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 895 observations removed because of NA values (LHS: 25, RHS: 895, Fixed-effects: 25).
> 
> tb2_2 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 901 observations removed because of NA values (LHS: 25, RHS: 901, Fixed-effects: 25).
> 
> tb2_3 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 915 observations removed because of NA values (LHS: 25, RHS: 915, Fixed-effects: 25).
> 
> tb2_4 <- feols(p3_CLI_DV ~  dem_lag1 * top1Share_interp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 +
+                  crisis_3_percent_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 916 observations removed because of NA values (LHS: 25, RHS: 916, Fixed-effects: 25).
> 
> 
> ## Gen. Table B2 ----------------------------------------------------------
> 
> coefNames <- c("dem_lag1" = "Democracy",
+                "top1Share_interp_lag1" = "Top 1\\% Income Share",
+                "dem_lag1:top1Share_interp_lag1" = "Top 1\\% Income Share $\\times$ Democracy",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list("(1)" = tb2_1, 
+                "(2)" = tb2_2, 
+                "(3)" = tb2_3, 
+                "(4)" = tb2_4)
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-b2.tex')
+ )
> 
> 
> # Figure B2 ---------------------------------------------------------------
> 
> # Uses estimates from Model 4, Table B2 above)
> cme_data <- plot_cme(tb2_4, 
+                      variables = "dem_lag1", 
+                      condition = "top1Share_interp_lag1",
+                      draw = FALSE)
> ggplot(cme_data,
+        aes(x = top1Share_interp_lag1, y = estimate)) +
+   geom_hline(yintercept = 0, color = 'red', size = .5, linetype = 'dashed') +
+   geom_line() +
+   geom_ribbon(aes(ymin = conf.low, 
+                   ymax = conf.high), alpha = .2) +
+   labs(x = 'Top 1% Income Share',
+        y = 'Marginal Effect of Democracy on CLI') +
+   theme_minimal(base_size = 14) +
+   theme(legend.position = 'none',
+         panel.grid.minor.x = element_blank(),
+         panel.grid.minor.y = element_blank())
> 
> ggsave(here('output', 'figures', 'figure-b2.pdf'),
+        width = 6, height = 4, units = 'in')
> 
> 
> # Table B3 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> ## Analyses ---------------------------------------------------------------
> 
> tb3_1 <- feols(p3_CLI_DV ~  dem_lag1 * gini_disp_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 1,009 observations removed because of NA values (LHS: 25, RHS: 1,009, Fixed-effects: 25).
> 
> tb3_2 <- feols(p3_CLI_DV ~  dem_lag1 * gini_disp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 1,011 observations removed because of NA values (LHS: 25, RHS: 1,011, Fixed-effects: 25).
> 
> tb3_3 <- feols(p3_CLI_DV ~  dem_lag1 * gini_disp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 1,024 observations removed because of NA values (LHS: 25, RHS: 1,024, Fixed-effects: 25).
> 
> tb3_4 <- feols(p3_CLI_DV ~  dem_lag1 * gini_disp_lag1 +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 +
+                  crisis_3_percent_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 1,024 observations removed because of NA values (LHS: 25, RHS: 1,024, Fixed-effects: 25).
> 
> 
> ## Generate Table B3 ------------------------------------------------------
> 
> coefNames <- c("dem_lag1" = "Democracy",
+                "gini_disp_lag1" = "Gini\\textsubscript{post-transfer}",
+                "dem_lag1:gini_disp_lag1" = "Gini\\textsubscript{post-transfer} $\\times$ Democracy",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list("(1)" = tb3_1, 
+                "(2)" = tb3_2, 
+                "(3)" = tb3_3, 
+                "(4)" = tb3_4)
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-b3.tex')
+ )
>
> # Table B4 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> setDT(df_t1)
> setorder(df_t1, code_cowc, p3)
> df_t1[ , polity_lag1 := shift(x = polity, n = 1L, fill = NA, type = 'lag'), by = code_cowc]
> 
> 
> ## Analyses ---------------------------------------------------------------
> 
> tb4_1 <- feols(p3_CLI_DV ~  polity_lag1 | code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 184 observations removed because of NA values (RHS: 184).
> 
> tb4_2 <- feols(p3_CLI_DV ~  polity_lag1 * top1Share_interp_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 869 observations removed because of NA values (RHS: 869).
> 
> tb4_3 <- feols(p3_CLI_DV ~  polity_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 875 observations removed because of NA values (RHS: 875).
> 
> tb4_4 <- feols(p3_CLI_DV ~  polity_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 890 observations removed because of NA values (RHS: 890).
> 
> tb4_5 <- feols(p3_CLI_DV ~  polity_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 +
+                  crisis_3_percent_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 891 observations removed because of NA values (RHS: 891).
> 
> 
> ## Gen. Table B4 ----------------------------------------------------------
> 
> coefNames <- c("polity_lag1" = "Polity",
+                "top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share",
+                "polity_lag1:top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share $\times$ Polity",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list("(1)" = tb4_1, 
+                "(2)" = tb4_2, 
+                "(3)" = tb4_3, 
+                "(4)" = tb4_4, 
+                "(5)" = tb4_5)
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-b4.tex')
+ )
> 
> # Table B5 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> ## Analyses ---------------------------------------------------------------
> 
> 
> tb5_1 <- feols(p3_CLI_DV ~  polyarchy_lag1 * top1Share_interp_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 838 observations removed because of NA values (LHS: 25, RHS: 838, Fixed-effects: 25).
> 
> tb5_2 <- feols(p3_CLI_DV ~  polyarchy_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 848 observations removed because of NA values (LHS: 25, RHS: 848, Fixed-effects: 25).
> 
> tb5_3 <- feols(p3_CLI_DV ~  polyarchy_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 885 observations removed because of NA values (LHS: 25, RHS: 885, Fixed-effects: 25).
> 
> tb5_4 <- feols(p3_CLI_DV ~  polyarchy_lag1 * top1Share_interp_lag1_ln +
+                  gdp_lag1_ln +
+                  gdppc_lag1_ln +
+                  trade_openness_lag1 +
+                  crisis_3_percent_lag1 | 
+                  code_cowc + p3,
+                data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 886 observations removed because of NA values (LHS: 25, RHS: 886, Fixed-effects: 25).
> 
> 
> ## Gen. Table B5 ----------------------------------------------------------
> 
> coefNames <- c("polyarchy_lag1" = "Polyarchy",
+                "top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share",
+                "polyarchy_lag1:top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share $\\times$ Polyarchy",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list(tb5_1, tb5_2, tb5_3, tb5_4)
> names(models) <- paste0("(", 1:length(models), ")")
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-b5.tex')
+ )
> 
> 
> # Table D1 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> ## Analyses ---------------------------------------------------------------
> td1_1 <- feols(p3_CLI_DV ~  dem_bmr_lag1 | code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 116 observations removed because of NA values (LHS: 25, RHS: 116, Fixed-effects: 25).
> 
> td1_2 <- feols(p3_CLI_DV ~  dem_bmr_lag1 * top1Share_interp_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 856 observations removed because of NA values (LHS: 25, RHS: 856, Fixed-effects: 25).
> 
> td1_3 <- feols(p3_CLI_DV ~  dem_bmr_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 862 observations removed because of NA values (LHS: 25, RHS: 862, Fixed-effects: 25).
> 
> td1_4 <- feols(p3_CLI_DV ~  dem_bmr_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 890 observations removed because of NA values (LHS: 25, RHS: 890, Fixed-effects: 25).
> 
> td1_5 <- feols(p3_CLI_DV ~  dem_bmr_lag1 * top1Share_interp_lag1_ln +
+               gdp_lag1_ln +
+               gdppc_lag1_ln +
+               trade_openness_lag1 +
+               crisis_3_percent_lag1 | 
+               code_cowc + p3,
+             data = df_t1[df_t1$p3 >= 21 & df_t1$top1_n >= 20, ])
NOTE: 891 observations removed because of NA values (LHS: 25, RHS: 891, Fixed-effects: 25).
> 
> ## Gen. Table D1 -----
> 
> coefNames <- c("dem_bmr_lag1" = "Democracy\\textsubscript{BMR}",
+                "top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share",
+                "dem_bmr_lag1:top1Share_interp_lag1_ln" = "ln Top 1\\% Income Share $\\times$ Democracy\\textsubscript{BMR}",
+                "gdp_lag1_ln" = "ln GDP",
+                "gdppc_lag1_ln" = "ln GDP per cap.",
+                "trade_openness_lag1" = "Trade Openness",
+                "crisis_3_percent_lag1" = "Economic Crisis"
+ )
> 
> models <- list(td1_1, td1_2, td1_3, td1_4, td1_5)
> names(models) <- paste0("(", 1:length(models), ")")
> 
> new_rows <- data.frame(matrix(c("Period Length", rep("3 years", length(models)),
+                                 "Country \\& Period FE?", rep("Yes", length(models)),
+                                 "Countries", getGroups(models)),
+                               ncol = length(models) + 1, byrow = TRUE))
> 
> attr(new_rows, 'position') <- c((length(coefNames) * 2) + 1:3)
> 
> msummary(models,
+          coef_map = coefNames,
+          stars = c('*' = .1, '**' = .05, '***' = .01),
+          add_rows = new_rows,
+          gof_map = list(list('raw' = 'nobs', 'clean' = 'Observations', 'fmt' = f),
+                         list('raw' = 'adj.r.squared', 'clean' = 'Adj. R$^2$', 'fmt' = '%.3f')),
+          notes = "Standard errors clustered by country. Independent variables lagged by 1 time period.",
+          escape = FALSE,
+          output = here('output', 'tables', 'table-d1.tex')
+ )
>
>
> # Figure D1 ---------------------------------------------------------------
> 
> ## Create subsets ---------------------------------------------------------
> 
> # Low inequality
> df_f3$lowInequality_subset_drop_bmr <- 1
> 
> for (r in which((df_f3$top1Share_interp >= median(df_f3$top1Share_interp, na.rm = T) |
+                  is.na(df_f3$top1Share_interp)) 
+                 & df_f3$dem_trans_bmr == 1)) {
+   df_f3[ r:(r+10), lowInequality_subset_drop_bmr := 0]
+ }
> 
> 
> # High inequality subsets
> 
> df_f3$highInequality_subset_drop_bmr <- 1
> 
> for (r in which((df_f3$top1Share_interp < median(df_f3$top1Share_interp, na.rm = T) |
+                  is.na(df_f3$top1Share_interp)) 
+                 & df_f3$dem_trans_bmr == 1)) {
+   df_f3[ r:(r+10), highInequality_subset_drop_bmr := 0]
+ }
> 
> data_bmr_low <- as.data.frame(df_f3[lowInequality_subset_drop_bmr == 1 & year >= 1945])
> data_bmr_high <- as.data.frame(df_f3[highInequality_subset_drop_bmr == 1 & year >= 1945])
> 
> pm_lowInequality_bmr <- PanelMatch(lag = 5, lead = 0:10, 
+                                    time.id = 'year', unit.id = 'cid', 
+                                    treatment = 'democracy_bmr', 
+                                    outcome.var = 'cli_overall_norm',
+                                    refinement.method = 'CBPS.weight',
+                                    covs.formula = ~ 
+                                      lag(cli_overall_norm_lm3, 1) +
+                                      lag(gdp_ln_lm3, 1) +
+                                      lag(gdppc_ln_lm3, 1) +
+                                      lag(trade_openness_lm3, 1) +
+                                      lag(top1Share_interp_ln_lm3, 1),
+                                    qoi = 'att', 
+                                    forbid.treatment.reversal = FALSE,
+                                    data = data_bmr_low)
Warning message:
In panel_match(lag, time.id, unit.id, treatment, refinement.method,  :
  non-numeric data exists. Only numeric (including binary) data can be used for refinement and calculations
> 
> pm_highInequality_bmr <- PanelMatch(lag = 5, lead = 0:10, 
+                                     time.id = 'year', unit.id = 'cid', 
+                                     treatment = 'democracy_bmr', 
+                                     outcome.var = 'cli_overall_norm',
+                                     refinement.method = 'CBPS.weight',
+                                     covs.formula = ~ 
+                                       lag(cli_overall_norm_lm3, 1) +
+                                       lag(gdp_ln_lm3, 1) +
+                                       lag(gdppc_ln_lm3, 1) +
+                                       lag(trade_openness_lm3, 1) +
+                                       lag(top1Share_interp_ln_lm3, 1),
+                                     qoi = 'att', 
+                                     forbid.treatment.reversal = FALSE,
+                                     data = data_bmr_high)
Warning message:
In panel_match(lag, time.id, unit.id, treatment, refinement.method,  :
  non-numeric data exists. Only numeric (including binary) data can be used for refinement and calculations
> 
> 
> ## Plot -------------------------------------------------------------------
> 
> nicePMplot(pm_lowInequality_bmr, 
+            data_bmr_low,
+            yupper = .4, ylower = -0.24,
+            seed = 1)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
$plot

$estimates
       estimate         2.5%     97.5%           5%        95%
t+0  0.02648984 -0.020680828 0.1147850 -0.017144844 0.09186636
t+1  0.07008537 -0.020122233 0.1804457 -0.011400169 0.15607936
t+2  0.05038417 -0.034277644 0.1587921 -0.027670420 0.13280147
t+3  0.04883744 -0.038928968 0.1631300 -0.030089316 0.13505910
t+4  0.10535175 -0.012271045 0.2315544  0.006548784 0.20753852
t+5  0.16013298  0.014082881 0.3348222  0.037444805 0.29425616
t+6  0.10538125 -0.010010462 0.2338711  0.002884703 0.20671987
t+7  0.12970527  0.001902489 0.2898149  0.015173577 0.25068878
t+8  0.11078584 -0.033564726 0.2932307 -0.019958575 0.24440600
t+9  0.09300072 -0.058185996 0.2744566 -0.043050774 0.22806136
t+10 0.08498934 -0.070279247 0.2698356 -0.054165432 0.21981764

> 
> ggsave(here('output', 'figures', 'figure-d1a.pdf'),
+        width = 6, height = 4, units = 'in')
> 
> nicePMplot(pm_highInequality_bmr, 
+            data_bmr_high,
+            yupper = .4, ylower = -0.24,
+            seed = 1)
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
Weighted Difference-in-Differences with Covariate Balancing Propensity Score
Matches created with 5 lags

Standard errors computed with 2000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
$plot

$estimates
        estimate        2.5%     97.5%          5%       95%
t+0  0.030729263 -0.05979052 0.1630097 -0.04253845 0.1308057
t+1  0.056962909 -0.09960841 0.2161153 -0.06175581 0.1829961
t+2  0.048503091 -0.11648040 0.2074633 -0.07467642 0.1732886
t+3  0.043489550 -0.12758183 0.2017236 -0.08322297 0.1689304
t+4  0.031269511 -0.15461512 0.1903008 -0.10098090 0.1557140
t+5  0.004211078 -0.20838408 0.1666562 -0.15107726 0.1357502
t+6  0.036042843 -0.19695390 0.2175509 -0.14151420 0.1877751
t+7  0.036001208 -0.21722059 0.2057667 -0.15065703 0.1741035
t+8  0.030498903 -0.22625000 0.1983880 -0.15597127 0.1681525
t+9  0.026140060 -0.23019415 0.1930420 -0.15966000 0.1633366
t+10 0.013444669 -0.24846073 0.1786200 -0.18293178 0.1515952

> 
> ggsave(here('output', 'figures', 'figure-d1b.pdf'),
+        width = 6, height = 4, units = 'in')
> 
> 
> # Figure E1 ----------------------------------------------------------------
> 
> # Done in stata. See ABC
> 
> # Figure F1 ----------------------------------------------------------------
> 
> ## Data ----
> df_t1 <- haven::read_dta(here('data', 'cli_p3.dta'))
> 
> ## Analyses ---------------------------------------------------------------
> 
> m <- interflex(estimator = 'binning',
+                Y = 'p3_CLI_DV',
+                D = 'dem_lag1',
+                X = 'top1Share_interp_lag1_ln',
+                FE = c('code_cowc', 'p3'),
+                Z = c('gdp_lag1_ln', 'gdppc_lag1_ln', 
+                      'trade_openness_lag1', 'crisis_3_percent_lag1'),
+                na.rm = TRUE,
+                data = as.data.frame(df_t1[df_t1$dem_lag1 %in% 0:1, ]))
Baseline group not specified; choose treat = 0 as the baseline group. 
> 
> plot(m, 
+      interval = c(1.5, 3.5),
+      Ylabel = "CLI",
+      Dlabel = "Democracy",
+      Xlabel = "ln Top 1% Income Share",
+      theme.bw = TRUE,
+      show.grid = FALSE,
+      cex.lab = .75,
+      cex.axis = .75)
> 
> ggsave(here('output', 'figures', 'figure-f1.pdf'),
+        width = 5, height = 4, units = 'in')
> 
> # Figure F2 ----------------------------------------------------------------
> 
> set.seed(94305)
> m_kern <- interflex(estimator = 'kernel',
+                     Y = 'p3_CLI_DV',
+                     D = 'dem_lag1',
+                     X = 'top1Share_interp_lag1_ln',
+                     FE = c('code_cowc', 'p3'),
+                     Z = c('gdp_lag1_ln', 'gdppc_lag1_ln', 
+                           'trade_openness_lag1', 'crisis_3_percent_lag1'),
+                     na.rm = TRUE,
+                     data = as.data.frame(df_t1[df_t1$dem_lag1 %in% 0:1, ]))
Baseline group not specified; choose treat = 0 as the baseline group. 
Cross-validating bandwidth ... 
Parallel computing with 4 cores...
Optimal bw=0.0957.
Number of evaluation points:50
Parallel computing with 4 cores...
> 
> plot(m_kern, 
+      interval = c(1.5, 3.5),
+      Ylabel = "CLI",
+      Dlabel = "Democracy",
+      Xlabel = "ln Top 1% Income Share",
+      theme.bw = TRUE,
+      show.grid = FALSE,
+      cex.lab = .75,
+      cex.axis = .75)
> 
> ggsave(here('output', 'figures', 'figure-f2a.pdf'),
+        width = 5, height = 4, units = 'in')
> 
> 
> # Figures G1 + G2 ------------------------------------------------------
> 
> ## Low Cov Bal ---- 
> low_refined <- as.data.table(get_covariate_balance(pm_lowInequality$att,
+                                                    covariates = c(
+                                                      'cli_overall_norm_lm3',
+                                                      'gdp_ln_lm3',
+                                                      'gdppc_ln_lm3',
+                                                      'trade_openness_lm3',
+                                                      'top1Share_interp_ln_lm3'),
+                                                    use.equal.weights = F,
+                                                    data= data_polity_low,
+                                                    plot = F))
> 
> low_unrefined <- as.data.table(get_covariate_balance(pm_lowInequality$att,
+                                                      covariates =  c(
+                                                        'cli_overall_norm_lm3',
+                                                        'gdp_ln_lm3',
+                                                        'gdppc_ln_lm3',
+                                                        'trade_openness_lm3',
+                                                        'top1Share_interp_ln_lm3'),
+                                                      use.equal.weights = T,
+                                                      data= data_polity_low,
+                                                      plot = F))
> 
> low_refined <- melt(low_refined,
+                     variable.name = 'variable',
+                     value.name = 'sd')
Warning message:
In melt.data.table(low_refined, variable.name = "variable", value.name = "sd") :
  id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns []. Consider providing at least one of 'id' or 'measure' vars in future.
> low_refined$t <- rep(-5:-1, 5)
> low_refined$refined <- 'Refined'
> 
> low_unrefined <- melt(low_unrefined,
+                       variable.name = 'variable',
+                       value.name = 'sd')
Warning message:
In melt.data.table(low_unrefined, variable.name = "variable", value.name = "sd") :
  id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns []. Consider providing at least one of 'id' or 'measure' vars in future.
> low_unrefined$t <- rep(-5:-1, 5)
> low_unrefined$refined <- 'Unrefined'
> 
> plotData_low <- rbind(low_refined, low_unrefined)
> 
> plotData_low$refined_p <- factor(plotData_low$refined, levels = c('Unrefined', 'Refined'))
> 
> ## High Cov Bal ----
> high_refined <- as.data.table(get_covariate_balance(pm_highInequality$att,
+                                                     covariates = c(
+                                                       'cli_overall_norm_lm3',
+                                                       'gdp_ln_lm3',
+                                                       'gdppc_ln_lm3',
+                                                       'trade_openness_lm3',
+                                                       'top1Share_interp_ln_lm3'),
+                                                     use.equal.weights = F,
+                                                     data= data_polity_high,
+                                                     plot = F))
> 
> high_unrefined <- as.data.table(get_covariate_balance(pm_highInequality$att,
+                                                       covariates = c(
+                                                         'cli_overall_norm_lm3',
+                                                         'gdp_ln_lm3',
+                                                         'gdppc_ln_lm3',
+                                                         'trade_openness_lm3',
+                                                         'top1Share_interp_ln_lm3'),
+                                                       use.equal.weights = T,
+                                                       data= data_polity_high,
+                                                       plot = F))
> 
> high_refined <- melt(high_refined,
+                      variable.name = 'variable',
+                      value.name = 'sd')
Warning message:
In melt.data.table(high_refined, variable.name = "variable", value.name = "sd") :
  id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns []. Consider providing at least one of 'id' or 'measure' vars in future.
> high_refined$t <- rep(-5:-1, 5)
> high_refined$refined <- 'Refined'
> 
> high_unrefined <- melt(high_unrefined,
+                        variable.name = 'variable',
+                        value.name = 'sd')
Warning message:
In melt.data.table(high_unrefined, variable.name = "variable", value.name = "sd") :
  id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns []. Consider providing at least one of 'id' or 'measure' vars in future.
> high_unrefined$t <- rep(-5:-1, 5)
> high_unrefined$refined <- 'Unrefined'
> 
> plotData_high <- rbind(high_refined, high_unrefined)
> 
> plotData_high$refined_p <- factor(plotData_high$refined, levels = c('Unrefined', 'Refined'))
> 
> ## Plot ---- 
> 
> ggplot(plotData_low, aes(y = sd, x = t, color = variable)) +
+   geom_line() +
+   geom_hline(yintercept = 0, linetype = 'dashed') +
+   scale_color_discrete(labels=c('CLI', 
+                                 'GDP',
+                                 'GDP per cap.',
+                                 'Trade Dep.',
+                                 'Top 1%')) +
+   scale_y_continuous(limits = c(-.5, 6.75),
+                      breaks = 0:7) +
+   facet_wrap(~refined_p) +
+   labs(y = 'SD',
+        x = 'Time to Democratization') +
+   theme_cowplot() +
+   theme(legend.position = 'bottom',
+         legend.title = element_blank())
> 
> ggsave(here('output', 'figures', 'figure-g1.pdf'),
+        height = 4, width = 7, units = 'in')
> 
> ggplot(plotData_high, aes(y = sd, x = t, color = variable)) +
+   geom_line() +
+   geom_hline(yintercept = 0, linetype = 'dashed') +
+   scale_color_discrete(labels=c('CLI', 
+                                 'GDP',
+                                 'GDP per cap.',
+                                 'Trade Dep.',
+                                 'Top 1%')) +
+   scale_y_continuous(limits = c(-.5, 6.75),
+                      breaks = 0:7) +
+   facet_wrap(~refined_p) +
+   labs(y = 'SD',
+        x = 'Time to Democratization') +
+   theme_cowplot() +
+   theme(legend.position = 'bottom',
+         legend.title = element_blank())
> 
> ggsave(here('output', 'figures', 'figure-g2.pdf'),
+        height = 4, width = 7, units = 'in')